Attention select de mass (dépendance de mixOmics) prend le pas sur dplyr
tCPMG_vector <- unlist(as.vector(t(CPMG)))
CPMG_vector <- unlist(as.vector(CPMG))
tNOESY_vector <- unlist(as.vector(t(NOESY)))
NOESY_vector <- unlist(as.vector(NOESY))
par(mfrow = c(2, 2))
hist(CPMG_vector, breaks = 100, xlab = 'CPMG', main = "Histogram of CPMG" )
hist(log2(CPMG_vector), breaks = 100, xlab = 'log2 CPMG', main = "Histogram of log2 CPMG")
hist(NOESY_vector, breaks = 100, xlab = 'Noesy', main = "Histogram of NOESY")
hist(log2(NOESY_vector), breaks = 100, xlab = 'log2 Noesy', main = "Histogram of log2 NOESY")df_row_sum <- apply(CPMG,MARGIN = 1, FUN= mean)
df_col_sum <- apply(CPMG,MARGIN = 2, FUN= mean)
meanr <- mean(df_row_sum)
meanc <- mean(df_col_sum)
minr <- min(df_row_sum)
maxr <- max(df_row_sum)
minc <- min(df_col_sum)
maxc <- max(df_col_sum)
medr <- median(df_row_sum)
medc <- median(df_col_sum)
df <- data.frame(
nb_ech = c(846,202),
moy = c(meanr,meanc),
median = c(medr,medc),
min = c(minr,minc),
max = c(maxr, maxc)
)
rownames(df) <- c("Echantillon","Bucket")
kable(df, align = "c", caption = "Tableau de moyenne des échantillons et buckets pour CPMG") %>%
kable_styling(
#bootstrap_options = "striped",
full_width = F,
position = "center")| nb_ech | moy | median | min | max | |
|---|---|---|---|---|---|
| Echantillon | 846 | 0.0059508 | 0.0059521 | 0.0058323 | 0.0060026 |
| Bucket | 202 | 0.0059508 | 0.0015731 | 0.0000588 | 0.1327695 |
tCPMG <- t(CPMG)
log2_tCPMG <- log2(tCPMG + 1e-8)
log2_CPMG <- log(CPMG + 1e-8)
###Standardisation par la médiane
stdr_medians <- apply(log2_tCPMG, 2, median)
series_median <- median(stdr_medians)
log2_CPMG_centered <- data.frame(matrix(
nrow = nrow(log2_tCPMG),
ncol = ncol(log2_tCPMG)))
colnames(log2_CPMG_centered) <- colnames(log2_tCPMG)
rownames(log2_CPMG_centered) <- rownames(log2_tCPMG)
for (j in 1:ncol(log2_tCPMG)) {
log2_CPMG_centered[, j] <- log2_tCPMG[, j] - stdr_medians[j] + series_median
}
### Mise à l'échelle par la plage interquartile. Si sd ne change rien
sampleIQR <- apply(log2_CPMG_centered, 2, IQR)
seriesIQR <- median(sampleIQR)
scalingFactors <- sampleIQR / seriesIQR
CPMG_standard <- data.frame(matrix(
nrow = nrow(log2_tCPMG),
ncol = ncol(log2_tCPMG)))
colnames(CPMG_standard) <- colnames(log2_CPMG_centered)
rownames(CPMG_standard) <- rownames(log2_CPMG_centered)
for (j in 1:ncol(log2_CPMG_centered)) {
CPMG_standard[, j] <-
(log2_tCPMG[, j] - stdr_medians[j] ) / scalingFactors[j] + series_median
}
par.ori <- par(no.readonly = TRUE)
par(mfrow = c(2, 2))
par(mar = c(5, 5, 2, 1))
#par(mar = c(5, 6, 4, 1)) ou 4 6 5 1
sample_size <- 30
## select sample indices
selected_samples <- sort(sample(
x = 1:ncol(tCPMG),
size = sample_size,
replace = FALSE))
boxplot(tCPMG[,selected_samples],
horizontal = TRUE,
yaxt="n",
las = 1,
main = "original values",
xlab = "value",
ylab = "sample")
boxplot(log2_tCPMG[,selected_samples],
horizontal = TRUE,
yaxt="n",
las = 1,
main = "log2-transformed",
xlab = "log2(value)",
ylab = "sample")
boxplot(log2_CPMG_centered[,selected_samples],
horizontal = TRUE,
yaxt="n",
las = 1,
main = "Median-based centered",
xlab = "log2(value)",
ylab = "sample")
boxplot(CPMG_standard[,selected_samples],
horizontal = TRUE,
las = 1,
yaxt="n",
main = "Standardized\n(median centering, IQR scaling)",
xlab = "log2(value)",
ylab = "sample")Quand transformations uniquement en log2 Les 2 premiers axes expliquaient plus de 62% de la variance exprimée.
Désormais ils n’expliquent plus que 42% de la variance exprimée.
CPMG_stdr <- t(CPMG_standard)
res.pca <- PCA(CPMG_stdr, graph = FALSE)
fviz_eig(res.pca, addlabels = TRUE, ylim = c(0, 40))Quelque soit le groupe étudié (sexe, âge, tabac), sur les axes 1 et 2, les différents points se superposent, il n’est donc pas possible de distinguer un groupe d’un autre.
L’individu L36R338 semble excentré sur l’axe3 (axe1-20) par rapport aux autres points et dans une moindre mesure V351R356 qui lui ne ressort pas sur les 2 premiers axes.
Contribution des 20 buckets qui contribuent le plus à la variance expliquée.
Cette contribution est >1
fviz_pca_biplot(res.pca, repel = TRUE, axes = c(1,2),
select.var = list(contrib = 20),
geom.ind = "point",
fill.ind = Sample_ID$SEX,
pointshape = 21, pointsize = 2,
palette = "jco",
addaEllipses = TRUE,
col.var = "contrib",
gradient.cols = "BuGn",
legend.title = list(fill = "Sexe", color = "Contrib"))Juste pour montrer le code mais les résultats sont identiques. Remarque : Chaque axe est symétrique donc qu’un individu soit à-20ici et +20 avec factominer revient au même
trans.pca <- pca(CPMG_stdr, ncomp=5, center = TRUE, scale = TRUE)
plotIndiv(trans.pca, comp = c(1,2),ind.names = FALSE,
group = as.numeric(as.factor(Sample_ID$SEX)),
legend.title = "sexe")Regarde si on peut observer naturellement une clusterisation en fonction d’une catégorie
Quelque soit la catégorie sexe, age tabac, nous n’observons aucune clusterisation.
res.spca.CPMG <- spca(CPMG_stdr, scale=TRUE,
ncomp=3, keepX=c(10,10,10))
plotIndiv(res.spca.CPMG, ind.names = FALSE,
group=Sample_ID$SEX,
pch = as.numeric(Sample_ID$SEX),
pch.levels =Sample_ID$SEX,
legend = TRUE)voir données NOESY
trans.plsda <- plsda(CPMG_stdr, Sample_ID$SEX ,ncomp=5, scale = TRUE)
plotIndiv(trans.plsda, ind.names = FALSE, legend=TRUE,
ellipse = TRUE, title="PLS-DA sur CPMG")Ne sont représentés uniquement les buckets qui contribuent à au moins 50%
plotVar(trans.plsda, var.names = FALSE, cutoff = 0.5,
title = "Correlation circle plot avec cutoff 0.5")Possibilité de faire de la prédiction dans le cas ou les ind seraient en sous groupes, ce qui n’est pas notre cas. ceci est juste un test
trans.testplsda <- plsda(CPMG_stdr, Sample_ID$SEX ,ncomp=5, scale = TRUE)
background <- background.predict(trans.testplsda, comp.predicted=2,
dist = "max.dist")
plotIndiv(trans.testplsda, ind.names = FALSE, legend=TRUE,
ellipse = TRUE, title="PLS-DA sur CPMG", background = background,
legend.title="sexe")Les variables sélectionnées pour la composante 1 sont exprimées de manière positives pour le sexe1 et négative pour le sexe2. Pour le statut tabagique la composante est exclusivement négatives pour le statut non fumeur.
MyResult.plsda2 <- plsda(CPMG_stdr, Sample_ID$SEX, ncomp=5)
plotLoadings(MyResult.plsda2, contrib = 'max', method = 'mean')Résultats simlaires à la PLSDA
trans.splsda <- splsda(CPMG_stdr, Sample_ID$SEX ,ncomp=5, scale = TRUE)
plotIndiv(trans.splsda, ind.names = FALSE, legend=TRUE, legend.title="sexe",
ellipse = TRUE, title="Sparse PLS-DA sur CPMG")trans.splsda <- splsda(CPMG_stdr, Sample_ID$AGE ,ncomp=5, scale = TRUE)
plotIndiv(trans.splsda, ind.names = FALSE, legend=TRUE, legend.title="age",
ellipse = TRUE, title="Sparse PLS-DA sur CPMG")trans.splsda <- splsda(CPMG_stdr, Sample_ID$tabac ,ncomp=5, scale = TRUE)
plotIndiv(trans.splsda, ind.names = FALSE,legend=TRUE, legend.title="Tabac",
ellipse = TRUE, title="Sparse PLS-DA sur CPMG")Représentation par sexe uniquement
Ne sont représentés uniquement les buckets qui contribuent à au moins 50%
trans.splsda <- splsda(CPMG_stdr, Sample_ID$SEX ,ncomp=5, scale = TRUE)
plotVar(trans.splsda, var.names = FALSE, cutoff = 0.5,
title = "Correlation circle plot avec cutoff 0.5")Représentation des 10 premières contributions
MyResult.splsda2 <- splsda(CPMG_stdr, Sample_ID$SEX, keepX=10)
plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean', legend.title= "Sexe")MyResult.splsda2 <- splsda(CPMG_stdr, Sample_ID$tabac, keepX=10)
plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean', legend.title= "Tabac")MyResult.splsda2 <- splsda(CPMG_stdr, Sample_ID$AGE, keepX=10)
plotLoadings(MyResult.splsda2, contrib = 'max', method = 'mean', legend.title= "Age")